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ABSTRACT 

The problem of the resolution of turbulent flows in adaptive mesh refinement (AMR) 
simulations is investigated by means of 3D hydrodynamical simulations in an idealised 
setup, representing a moving subcluster during a merger event. AMR simulations per- 
formed with the usual refinement criteria based on local gradients of selected vari- 
ables do not properly resolve the production of turbulence downstream of the cluster. 
Therefore we apply novel AMR criteria which are optimised to follow the evolution of 
a turbulent flow. We demonstrate that these criteria provide a better resolution of the 
flow past the subcluster, allowing us to follow the onset of the shear instability, the 
evolution of the turbulent wake and the subsequent back-reaction on the subcluster 
core morphology. We discuss some implications for the modelling of cluster cold fronts. 

Key v^rords: Hydrodynamics ~ Instabilities - Methods: numerical - Galaxies: clusters: 
general - Turbulence 



1 INTRODUCTION 

The role of turbulence in astrophysics and, in particu- 
lar, in cosmic structure formation has been attracting in- 
creasing attention, along with observational and theoreti- 
cal indications of its relevance. There are different obser- 
vational hints converging towards t he turbulent character 
of the intra-cluster medium (ICM) (IS chuec ker et al.l 20041: 
Churazov et~aLl |2004 lEnfilin fc Vogf 2006; Rebusco et all 



20051 . l2006l . but see also lFabian et aLlbOOa l . From an obser- 



vational point of view, the next generation of X-ray observa- 
tories will be able to observe the Doppler shifts of emission 
lines of heavy ions due to bulk motions in the ICM, thus 
providing information about the turbulent state of the flow. 
Numerical studies predict that gas bulk motions and tur- 
bule nce in the ICM can be a consequence of cluster merg- 
ing llRickerlll994 iNorman fc Brvaiil Il999l : iTakizawal |2000| ; 



Ricker fc Sarazinll200ll: iDolag et aLlbOOSl ). In their model 



Subramanian et al.l ( 20061 ) identify three physical regimes 



for turbulence production and decay in clusters, the latest 
one being dominated by turbulent production in the wakes 
of minor mergers. This phase would play a key role also for 
the magnetic fleld amplification in the ICM. 

Merging is a crucial process in cosmic structure for- 
mation. It is predicted in the framework of hierarchi- 
cal clustering and can be stu died theoretically by means 
of numerical simulati ons (e.g. iRoettiger et al. r il996l . 1 19981 : 
[McCarthy et al.ll2007h . Observationally, modern X-ray ob- 
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servatories are able to detec t substructures in the ICM 
that are indicative of merging l|Sunvaev et al.|[2003t ). Among 
the first results of galaxy cluster observations with Chan- 
dra was the discovery of "cold fronts" in the merging clus- 
ters A2142 and A3667, and later in many other clusters 
(see iMarkevitch fc VikhlininI |2007| for a review). The ori- 
gin of cold fronts in merger clusters has been ascribed 
to infalling substructures, although the cold fronts in the 
cores of cooling core clusters have a different interpretatio n 
(|Tittlev fc Henriksenllioosi : lAscasibar fc Markevitchll2006l ). 

Cosmological simulations of galaxy cluster formation 
have confirmed the merging scenario for the transient for- 
matio n of cold fronts (iBialek et al.ll20o3 : iNagai fc Kravtsovl 
l2003h . The problem has also been addressed with sim- 
plified setups, which allow a better control over the in- 

jAcreman et al.j |2003l : 
. al.. ,20071) and 



iHeinz et all 12003 


: Asai et al.l 20041: 


3D simulations ( 


Takizawal l2005allbl: 


iDursi fc Pfrommei 20081). 
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The theoretical study of turbulence in the framework 
of full cosmological simulations is a challenging task and 
presents the typical problems of numerical simulations of 
strongly clumped media. Adaptive mesh refinement (AMR) 
is a viable tool for saving computational resources and mod- 
elling the large dynamic range in a proper way (|Normanl 
l2005h . The choice of the most suitable mesh refinement cri- 
terion is thus a delicate compromise between following the 
fiow structure in the most accurate way and exploiting the 
advantage of saving computational time and memory. 

This study (Paper I) is the first part of a broader project 
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on the physics of galaxy clusters, focused on investigat- 
ing the generation of turbulence in the ICM by means of 
AMR simulations. In Paper I, we present new AMR crite- 
ria that are more suitable for resolving the production of 
turbulence than those commonly used in cosmological simu- 
lations. They were tested in a simple setup, reminiscent of a 
wind tube experiment, but designed in order to be represen- 
tative of an idealised subcluster merger. This test problem 
is a useful and controlled testbed for the new tools, far from 
the complexity of a more realistic cosmological simulation. 
Nevertheless, it is similar to previous studies of subcluster 
mergers and can provide useful insights about the physics of 
the minor mergers, relevant for the generation of turbulence 
in the ICM. In contrast with earlier work, we focus our anal- 
ysis on the turbulent wake of the subcluster rather than on 
the cold front, showing the difference of the evolution of the 
Kelvin-Helmholtz instability (KHI) resulting from the use 

of different AMR criteria. 

The second step, described in llapichino fc Niemeveij 

l|2008 l) (hereafter Paper II), is to apply the same tools to 
simulations of structure formation and to study the turbu- 
lent flow in the ICM. The results of the present work will 
thus be compared with the outcomes of more realistic sim- 
ulations of minor mergers. 

The paper is structured as follows: the new refinement 
criteria are defined in in Sec. (2] The numerical setup of the 
performed simulations and the features of the KHI are in- 
troduced in Sects. [3] and m respectively. The results are pre- 
sented in Sec. [5] and discussed in Sec. [6] 



2 REFINEMENT CRITERIA AND 

RESOLUTION OF TURBULENT FLOWS 

Three different AMR criteria are used in the simulations 
discussed in this work. A rather general and widely usecj^ 
criterion is based on the local gradients of all thermodynam- 
ical variables. If the relative slope \{q{i-\-l) — q{i-'l)) / q{i)\ of 
a variable g in a cell i is larger than a given threshold, that 
cell is marked for refinement. Numerical tests showed that, 
for our subcluster problem, this criterion quickly leads to a 
high level of refinement in the entire computational domain. 
This shortcoming was fixed by allowing the refinement on 
the local gradients of only selected variables. The AMR cri- 
terion named "1" in our work is based on the local gradients 
of density and internal energy, with thresholds 0.24 and 0.25 
respectively, tuned for optimal performance. 

Two novel AMR crite ria, developed for trac king the evo- 
lution of a turbulent flow (|Schmidt et al.|[2008h . were tested 
in this study. In both cases, the control variables for re- 
finement are given by scalars probing small-scale features 
of the flow. One example is the modulus of the vorticity 
vector ijj = SI X V (the curl of the velocity field) that is 
expected to become high in regions filled by turbulent ed- 
dies. In addition, the mechanism for triggering refinement 
has been modified. Rather than normalising the control vari- 
ables in terms of characteristic quantities and comparing to 
prescribed threshold values, the new criteria use regional 

^ This is the criterion used in public code ENZO, employed for 
this work (Sec. |3]l. 



Table 1. Scheme of the performed simulations. 



Simulation 


number of AMR levels 


AMR criterion 


A 


5 


1 


B 


4 


1 


C 


6 


1 


D 


5 


2 


E 


5 


2+3 



thresholds for triggering refinement which are based on a 
comparison of the cell value of the variable q{x,t) with the 
average and the standard deviation of q, calculated on a 
local grid patch: 

q{x,t)^ {q)^{t)+aX^{t) (1) 

where Xi is the maximum between the average (q) and the 
standard deviation of q in the grid patch i, and a is a tun- 
able parameter. This technique can easily handle highly in- 
homogeneous problems such as subcluster mergers without 
a priori knowledge of the flow properties. 

We define criterion "2" using the square of the vortic- 
ity, (jj^, as the control variable. Refining upon criterion "2" 
only, one can expect to resolve turbulent eddies which are 
associated with high vorticity but not flow features arising 
from pure compression such as shock fronts. For this reason, 
in criterion "3" , the control variable is given by the rate of 
compression, i.e., the negative time derivative of the diver- 
gence d = V-u. In combination with refinement by vorticity, 
this criterion is intended to refine flow regions in which steep 
density gradients are developing. One could also refine by —d 
or in order to capture compression effects. However, this 
would give rise to full refinement of highly compressed gas, 
despite the fact that the interiors of compressed regions do 
not necessarily exhibit rich small-scale structure. Maxima 
of the rate of compression, on the other hand, are typically 
found at the interfaces between dense and rarefied gas. As we 
have seen in applications, a drawback of criterion "3" is that 
refinement can be triggered prematurely in some instances. 
This undesired side-effect can be suppressed by introducing 
a lower cutoff for the refinement thresholds. 



3 NUMERICAL SIMULATIONS 

The list of the performed calculations is reported in Table 
[1] In the first column, the simulations are identified by a 
letter. The second column reports the number n of AMR 
levels which are allowed in the simulation. The root grid 
resolution of the runs is only 16"^ but the effective resolution 
is finer because AMR is used. The refinement factor is 2, 
thus the effective resolution is (16 x 2")^ grid cells. In the 
standard case of five additional AMR levels, the effective 
resolution is therefore (16 x 2^)^ = 512^ grid elements. 

The third column of Table [T] indicates which criterion 
was used for performing the grid refinement, according to the 
definitions of Sec. [5] The simulations A, B and C make use of 
criterion "1" in a resolution study, whereas the simulations 
D and E were performed with the novel AMR criteria "2" 
with threshold a2 = 0.2, and "2-|-3" with 02 = 1.0, 03 = 2.0, 
respectively. 

The 3D simulations were performed using the AMR, 
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:rid-based hybrid (h ydrodynamical plus N-Body) code enzo 



(o'Shea et"ai1l2005l n. with the N-Body section of the code 
switched off. The primary hydrodynamic method used in 
ENZO is based on the piece wise parabolic method (PPM) 
I Woodward fc Colelia 1984 modified for the study of 
cosmology IjBrvan et al. I Il995l ). Also available is the hy- 
drodynamic fi nite-difference algorithm of the ZEUS code 
ijStone &i Norm an 1992a,j3), used for tests in Sec. I5.3.T1 

Since we assumed ideal adiabatic physics for our cal- 
culations, the simulations are scale-free and all quantities 
were rescaled to get nu mbers of the order unity, as shown in 
I Robinson et al.l l|2004l ). However, for the sake of clarity, all 
quantities are reported in CGS units. An adiabatic equation 
of state was used, with 7 = 5/3. 

In the initial conditions, an isothermal and spherical 
symmetric object is set in the co mputational domain accord- 
ing t o a beta density profile (ICavaliere fc Fusco-Femiand 



p{r) = pc 



1 + 



-3/3/2 



(2) 



where p is the gas density, r is the distance from the cen- 
tre, pc is the central density, and Vc is the core radius. The 
most important difference between this setup a nd a typical 
wind tube experiment (cf. iMurrav et al. l ll993l : IVietrr et al.l 
ll997l : [A"gertz et al]l2007l ) is the presence of gravity. The grav- 
itational acceleration is provided by a fixed spherical dark 
matter distribution following a King profile. 



Pdm(j') = pdm.c 



1 + 



-3/2 



(3) 



with Pdm,c = 10 Pc- A cutoff was set at r = 5 r^, both for 
the dark m atter density pro file and the gravitational accel- 
eration (cf. lTakizawall2005al ) . Given the King profile for the 
dark matter component of an isothermal cluster, the beta 
profile is the hydrostatic equilibrium solution for the gas 
density. The hydrostatic equilibrium is enforced on the dis- 
crete computational grid by the relaxat i on pr ocedure of the 
initial state described bv lZingale et al.l (|2002h . 

The interaction of the subcluster with a background 
wind mimics the merger process in a simplified manner. 
A uniform velocity field along the x-axis is imposed in the 
background medium. In this simplified setup, we implicitly 
assume that the merging process has a negligible impact on 
the dark matter density profile of the subcluster. 

The numerical values of the parameters defined above 
were chosen as typical of a subcluster merger. We set pc = 
6.3 X 10"^^ g cm"^, rc = 7.7 x 10^^ cm = 250 kpc and 
P — 0.6, which result in a subcluster temperature kT = 
3.65 keV. The background medium has constant values of 
density pb ~ 7.9 x 10~^^ g cm~"^ and temperature fcTb = 
8.0 keV. The background velocity is chosen as Hb = 1.6 x 
10^ km s~^, Db/cb = 1.1, where Cb is the sound speed of 
the background. The interface between the subcluster and 
the uniform medium is defined at the locations where the 
subcluster and background gas pressures are equal. 

The computational domain has a size of [0,4 Mpc]"^, 
meaning that the effective resolution of 512"^ corresponds 



^ ENZO homepage; | http://lca.ucsd.edu/portaI/software/enzo] 
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to an effective spatial resolution of 7.8 kpc. The boundary 
conditions are inflowing at the left x-boundary and outflow- 
ing elsewhere. The subcluster is initially resolved with the 
finest AMR level available, in order to minimise any map- 
ping errors related to the hydrostatic equilibrium of the ini- 
tial state. 



4 THE KELVIN-HELMHOLTZ INSTABILITY 

The simulations are carried out in the rest frame of the 
subcluster, set at rest in a moving background medium. 
The outer subcluster interface is thus subject to the Kelvin- 
Helmholtz instability. If gravity is neglected, the KHI may 
grow for all values of wavenumber k. The linear growth 
timescale of the KHI, in the case of a flat interface be- 
tween two incompre ssib le fluids of differe nt densities, is 
(|Drazin fc Reidll2004l ; cf. lAgertz et~alll2007l ): 



Tkh = 



27r(pi + P2) 
k{pip2y^^V 



(4) 



The geometry of the subcluster is different from this ide- 
alised case, but an order-of-magnitude estimate of tkh can 
be obtained by setting pi = 0.3 pc (subcluster density at the 
interface with the background medium), p2 ~ pb, k = 2n/R, 
where i? 1.7 Tc is the distance from the subcluster centre 
to the interface with the background medium. The velocity 
V is smaller than Vh because a bow shock forms in front of 
the subcluster and V has to be d etermined from the post- 
shock conditions (cf. lAgertz et al.|[2007, ). For this rough esti- 
mate it is adequate to assume V ~ 0.7 Ub (also in agreement 
with the simulations). With these choices, tkh ~ 0.8 Gyr. 
This is an estimate of the time required to the KHI to per- 
turb the spherical shape of the subcluster substantially. We 
set the total simulation time to 4 Gyr. 

The previous timescale analysis motivates the choice 
of neglecting radiative cooling in the performed simulations. 
The effect of cooling is not expected to be prominent because 
both for the subcluster and the background medium the 
cooling timescale is about two times longer than the total 
simulation time, and much longer than tkh. 

The presence of the gravitational acceleratio n sta- 
bil ises a perturbation of wavenumber k if (|Patersonlll984l : 
cf. lMurrav et al.lll993D 



PiP2V''fc 



Pi 



(5) 



where g is the modulus of the gravitational acceleration at 
the interface. In the subcluster case, this stabilisation only 
affects wave numbers that are much smaller than 2tt/R. 
Gravity thus fails to stabilise the subcluster against the KHI. 

The loss of gas from the subcluster atmosphere due 
to the ram pressu re of the background flow, discussed by 
iHeinz et af] l|2003h . is not effective for our choice of param- 
eters. The gas stripping and its subsequent mixing into the 
ambient medium is instead caused by the KHI, and will be 
described in Sect. 15.31 
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t - 1 Gyr 



t - 2 Gyr 




t - 3 Gyr 



t = 4 Gyr 



Figure 1. Density slices (xz plane) of a part of the computational domain (2.6 X 1.6 Mpc), showing the subclustcr evolution at different 
times, for the run A. The density is linearly colour coded, following the colour bar on the upper left. Time is indicated at the lower left 
of each panel. 



5 RESULTS 



5.1 Evolution of the subcluster 

The morphological evolution of the idealised subcluster is 
shown in the four density slices of Fig. [1] The figure refers 
to run A, which (as discussed in Sec. 15. 2|) has a suitable 
spatial resolution for our problem. 

At t = 1 Gyr, a bow shock has formed in front of the 
subcluster and the KHI is growing at the sides. At t = 2 Gyr 
a turbulent, eddy- like flow in the subcluster wake is clearly 
visible. Yet it is not resolved at later times, as will be dis- 
cussed in Sec. 15.31 The KHI leads to mass stripping in the 
subcluster, which will be quantified in Sect. 15.21 

From t — 2 Gyr, the morphology which is visible 
at the left interface of the subcluster (Fig. [SJ closely re- 
sembles a typical cold fro nt structure. Based on Chan- 
dra observations of A3667, Vikhlinin et al.l l|2001al lbh and 
IVikhhnin fc Markevitclil (|200a 7 state that the very existence 
of the cold fronts requires the suppression of the KHI and of 
the thermal conduction at the interface between the sub- 
cluster gas and the ICM, likely cause d by the presenc e 
magnetic f i elds. T he MHD simulations of lAsai et al. 
lAsai et al.1 l|2005l ) and lAsai et al] (|2007| ) include anisotropic 
thermal conduction and show an effective stabilisation of 
the cold front. The suppression of the thermal conduction is 
needed to reproduce the observed temperature profiles of the 




Figure 2. Temperature slice {xz plane) of the subcluster a.t t = 
2 Gyr for the run A (AMR on gradients of density and internal 
energy). The temperature is linearly colour coded in gray-scale. 



cold f ront in A3667 in the SPH simulations of IXiang et al.l 
|2003). 

Our simulations incl ude neither heat c onduction nor 
magnetic fields. Similar to lHeinz et al.l (|2003l ) . a "cold front- 
like" structure is thus visible but, according to the above 
cited simulations, it does not catch the physics of the astro- 
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12 3 4 

time [GyrJ 



Figure 4. Temporal evolution of subcluster baryonic mass frac- 
tion for the runs A (solid line), B (dotted line) and C (dashed 
line). The mass fraction is normalised to its value at t = Gyr. 




Figure 5. Density slice as in Fig.[T] upper right panel, with the 
mesh structure superimposed. Grids of AMR levels from to 3 
are rendered as mesh structures, whereas for ease of visualisation 
grids of level 4 and 5 are only rendered with colours green and 
blue, respectively. 



physical problem. However, its analysis is not the primary 
goal of this work. In Sec. (6] we will extend the discussion of 
the subcluster evolution to examples which take into account 
such additional physics. 



5.2 Resolution study 

A comparison of AMR runs with different resolution is use- 
ful for understanding how effectively (or not) the idealised 
subcluster is modelled in the simulations. With this intent, 
besides the presented run A we performed the additional 
runs B and C (cf. Table [!}, which implement the same AMR 
criterion as A and have the same root grid resolution, but 
allow different numbers of additional AMR levels. Figure [3] 
permits a morphological comparison at t = 2 Gyr. The run 
B allows four additional levels of refinement and has a lower 
effective resolution than run A. Comparing Fig. O left, and 
the upper right panel of Fig. [1] a less developed shear in- 
stability and an almost regular pattern in the flow at the 
subcluster tail is noticeable in the former. The crucial fea- 
tures of the subcluster evolution are therefore not properly 
followed at this level of resolution. On the other hand, the 
analysis of more resolved run C (six additional AMR lev- 
els) shows obviously a better resolution than run A and an 
enhanced development of the KHI at the sides of the sub- 
cluster. Nevertheless, it does not provide any new element 
for the study of the problem, and therefore the effective res- 
olution of the run A is considered to be adequate for our 
investigation. 

A more quantitative comparison is obtained by defin- 
ing a diagnostic of gas stripping resulting from the develop- 
ment of the KHI. A suitable qu antity is adapted fr om the 
"cloud mass fraction" defined bv lAgertz et~aLl l|2007l ). hence 
we consider the gas with T < 0.9 Th and p > 0.32 pc as 
belonging to the subcluster. After a transient phase related 
to the development of the bow shock, where the "subclus- 
ter mass" seems to inc rease, this quantity ( Fig. |4]) decreases 
with time. Similar to lAgertz et al.l l|2007l ), the subcluster 
stripping is more pronounced in the more resolved runs, due 
to the better resolved small scale instability and mixing. 



I 




Figure 6. Slice of the subcluster at t = 2 Gyr for the run A 
(cf. Fig.[T] upper right panel), showing the square of the vorticity 
modulus uj'^. The quantity is colour coded on a logarithmic scale 
in code units. 



5.3 AMR study 

One of the most interesting features of the subcluster evo- 
lution in run A is the loss of resolution in the wake at late 
times. This is an effect of the AMR criterion chosen in that 
simulation, which cannot effectively refine the turbulent tail 
(or, at least, without wasting the advantage of AMR, with 
unnecessary refinement in almost the entire computational 
domain) . 

In a study of cold front physics, one could think that 
the refinement of run A (Fig. [SJ is acceptable, because the 
attention is focused on the better resolved part of the compu- 
tational domain. We dispute this statement in this section. 
To this aim, we checked whether alternative AMR criteria 
can better track the turbulent wake, provided the rest of the 
cluster morphology is also appropriately refined. 

Velocity fiuctuations at all scales are the prominent fea- 
ture of a turbulent flow. It suggests that quantities related to 
the spatial derivatives of velocity can be particularly suitable 
for the characterisation of the flow. Figure |6l which shows a 
slice of the square of the vorticity modulus in run A, clearly 
confirms this idea: the morphology of the turbulent fiow is 
well tracked in comparison with the less contrasted density 
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Figure 3. Density slices of the subcluster at t = 2 Gyr (cf. Fig. [T] upper right panel), in simulations with different effective resolution. 
Left: run B. Right: run C. The density is linearly colour coded, according to the colour bar on the left. 






Figure 7. Slices (xz plane) for the run D, showing density (left) and square of the vorticity modulus, with the mesh superimposed and 
visualised according to the rendering used in Fig. |5] (right). Upper panels: slices i = 2 Gyr. Lower panels: slices t = 3 Gyr. 



and temperature eddies, which are not very suitable for trig- 
gering the mesh refinement. 

The refinement criteria "2" and "3", presented in 
Sect. [21 are based on this concept. Their effect on the res- 
olution of the turbulent subcluster wake can be evaluated 
by the analysis of the runs D and E (Figs. [7]and[Sl respec- 
tively, at t = 2 and 3 Gyr). From the density and vorticity 
slices one can immediately recognise that, in both simula- 
tions, the subcluster wake is effectively resolved down to the 
finest available AMR levels. Moving out of the subcluster, 
the refinement level decreases gradually. 

A closer look reveals that there are some differences 



between the two simulations. As anticipated in Sec. (2] the 
refinement criterion "2" alone is unable to resolve the bow 
shock in front of the subcluster at the finest levels of refine- 
ment (run D). It is difficult to state whether this has any 
effect on the development of the KHI, but it is certainly not 
desirable to underresolve such an important feature of the 
problem, well resolved in run E. 

The subcluster tail in run £ at t = 2 Gyr (Fig. [S] upper 
left) has a more filamentary pattern, intriguingly similar to 
the more resolved run C (cf. Fig. O right). 

At t = 3 Gyr, in both runs D and E the subcluster 
appears more perturbed and prone to the KHI than run 
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Figure 8. Same as Fig. [7] for run E. 



A (Fig. [T] lower left panel). Since the subcluster front is 
well resolved with the refinement criteria used in all simu- 
lations, this difference is (at least partly) to be ascribed to 
the back-reaction of the tail. The turbulent eddies in the 
subcluster wake are better resolved in runs D and E, and 
partly disturb the morphology of the subcluster core even 
though it is well resolved. This effect of back-flow is rather 
similar to that described bv lHeinz et al. II2OO3), which iden- 
tify the displacement of the subcluster core with respect to 
the potential well as an additional source. 

As a measure of the action of the KHI, the subcluster 
mass has been introduced in Sec. 15.21 It is useful to inves- 
tigate quantitatively the effect of the back-flow, comparing 
the time evolution of this quantity in the runs A, C, D, 
and E (Fig|9|. Interestingly, after t = 2 Gyr the stripping 
for the runs D and E is more effective than for the run A, 
and follows more closely the curve of the more resolved run 
C. The seeming analogy comes from the convergence of two 
different features of the small-scale mixing, which for run C 
is more active at the sides of the subcluster, whereas for the 
runs D and E is triggered mainly by the tail back-flow. 

In order to further quantify the effectiveness of the 
tested AMR criteria, the mass-weighted root mean squared 
(rms) velocity iirms is calculated in a test volume of 
(400 kpc)^, located immediately downstream of the subclus- 
ter, aX t = 2) Gyr (Tabled - The simulations using the novel 
AMR criteria are able to get much higher turbulent veloci- 
ties than r un A. These values are in the range predicted the- 
oretically (jSubramanian et al. 20061) and found numerica lly 
in similar works ( Acreman et al.l 20031 : iTakizawal l2005al lbh . 
Interestingly, the new AMR criteria bring Vrms close to the 



s 0.55 



y. 1 ' ' ' ' 


1 11 11 1 11 11 _ 




run A 




run C 




run D - 




run E 


1 , , , , 





time [GyrJ 

Figure 9. Detail of the temporal evolution of normalised sub- 
cluster mass fraction for the runs A (solid line), C (dotted line), 
D (dashed line) and E (dot-dashed line). 



value which is obtained by an equivalent static grid simula- 
tion (cf. Sect.[5X2j. 

A closer analysis of the subcluster tail shows that, at a 
distance from the subcluster core r « 5 rc, the density struc- 
ture slightly loses contrast. This is an artifact introduced by 
the cutoff in the dark matter density and gravitational accel- 
eration (Sect.[3|. We tested that the results are not harmed 
by this shortcoming. Moreover, the vorticity pattern is basi- 
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Table 2. Mass- weighted rms velocity in the test volume of 
(400 kpc)^, at t = 3 Gyr. The run labelled as "Static" will be 
introduced in Sect. I5.3T21 



run 


t^rms (km S ^) 


A 


260 


D 


530 


E 


530 


Static 


620 



Table 3. Occupation fraction of the AMR levels at t = 2 Gyr, for 
the runs A, D and E (the AMR level corresponds to the root 
grid). The fraction is normalised to the whole computational do- 
main. For computational efficiency, in all simulations the refine- 
ment is allowed only in the region of interest around the subclus- 
ter, with a volume fraction 0.608. 



AMR level 


run A 


run D 


run E 


1 


0.150 


0.156 


0.422 


2 


0.100 


0.076 


0.255 


3 


0.052 


0.045 


0.126 


4 


0.027 


0.031 


0.074 


5 


0.015 


0.022 


0.045 



cally not affected by the cutoff (Figs. [7] and [S]) and the same 
is true for the new refinement criteria. 

From the point of view of the AMR performance, Table 
[3] summarises the volume occupation fractions at different 
levels for three of the runs under examination, ai t = 2 Gyr. 
As expected, the simulations D and E have larger occupa- 
tion fractions than run A. In all cases, the ratio between the 
occupation fractions at levels n and n+1 lies between 1.4 and 
2.0, and only a small fraction of the computational domain 
is refined at the higher levels. Though the overall quality of 
the AMR operation seems thus acceptable, a visual inspec- 
tion of Fig. [7]and[Sl right panels, shows that a more efficient 
use of refinement is made in run D than in run E. The latter 
refines the front bow shock better than the former but, de- 
spite of the implemented refinement cutoff, it triggers some 
spurious grids at the sides of the subcluster. This issue is 
not severe and is perfectly manageable within the available 
computational resources. In a general sense, this is a good 
example of tuning AMR to specific problems, finding a dif- 
ficult equilibrium between a accurate description of the fiow 
and a convenient use of the tool. 

5.3.1 Comparison with results from the ZEUS solver 

It is well known that the simulation of the development of 
hydrodynamical instabilities and the onset of turbulence de- 
pends sensitively on the adopte d numerical s c heme . This 
issue is extensively addressed by lAgertz et af] (|2007l ). m a 
simulation setup which differs from ours only for the role 
of gravity, and for the density profile of the subcluster. A 
similar, detailed study of our setup is out of the scope of 
this paper, but we performed a simulation similar to run A 
(5 AMR levels, AMR criterion "1" ) using th e zeus solver 
available in enzo (cf. IStone fc Normanlll992al lb[) instead of 
PPM. The u s e of e nzo with this solver has been tested by 
lAgertz et all (12007 1) in static grid simulations, whereas here 
it is applied to an AMR run. 




Figure 10. Density slice of the subcluster at t = 2 Gyr (cf. Fig.[T] 
upper right panel), in the simulation with the AMR parameters 
as in run A, but using the ZEUS solver of ENZO. The density is 
linearly colour coded following the colour bar on the left. 




12 3 

time [GyrJ 



Figure 11. Temporal evolution of subcluster baryonic mass frac- 
tion for the runs A (solid line) and the run using the ZEUS solver 
(dotted line). The mass fraction is normalised to its value at 
t = Gyr. 



The results of this simulation show clearly that the ZEUS 
solver is unable to follow the onset of the KHI and the sub- 
sequent development of turbulent motions in the subcluster 
wake, where no clear turnover eddy is visible (Fig. I10|l . The 
evolution of the subcluster mass fraction (Fig. Ilip is typical 
of simulations where the turbulent mixi ng is not efficiently 
acting (examples are in some figures of lAgertz et al.|[2007l . 
or in our run B). Moreover, both the magnitude and the 
turbulent features of the flow are not well reproduced in the 
wake immediately behind the substructure (compared with 
run A in Fig. [T2|) . 

The described features do not depend on the chosen 
artiflcial viscos ity of the ZEUS sol ver (set at 2.0 in the per- 
formed run; cf . lAgertz et allbOOTh . We compared the results 
with a static run performed with the ZEUS solver and the 
same level of resolution of the AMR run, finding no sub- 
stantial difference. The relevance of this test confirms that 
any conclusion on simulations of turbulent flows must be 
critically considered in light of the used numerical scheme, 
as will be further discussed in Sect. [6] 
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Figure 12. Slices of the subclustcr at t = 2 Gyr (cf. Fig.[T] upper right panel), showing the velocity magnitude. Left: run A. Right: run 
using the ZEUS solver. The velocity magnitude is linearly coded such that that the colour-bar ranges from v = (black) to the imposed 
background velocity (v = 1.6 X 10"^ km s~^, white). 




12 3 4 

time [GyrJ 



Figure 14. Temporal evolution of subcluster baryonic mass frac- 
tion for the runs A (solid line), the run with a 512^ static grid 
(dotted line) and the run D as representative of the runs with 
new AMR criteria (dashed line). The mass fraction is normalised 
to its value at t = Gyr. 



5.3.2 Comparison with a static grid simulation 

The computational benefit from the use of AMR is espe- 
cially needed in simulations where the corresponding static 
grid is too expensive; this is often the case in cosmologi- 
cal simulations of structure formation. On the other hand, 
most of the presented runs have an effective resolution of 
512'^, making the comparison with a static grid simulation 
still possible. In view of further and computationally more 
demanding applications, a comparison of AMR and static 
grid simulations is therefore particularly useful. 

A static grid simulation of the subcluster setup has been 
performed with a fixed resolution of 512^ cells, correspond- 
ing to the effective resolution of the runs A, D and E. The 
density slice at t = 2 Gyr (Fig. 1131 left) is similar to the 
AMR runs, with only a higher level of symmetric structure 
in the wake. The symmetry is not preserved at t = 3 Gyr 
(Fig. 1131 right), when the turbulent tail is very similar to 
the runs with new AMR criteria (Figs. [7]and[8)|. 

The evolution of the subcluster mass fraction, shown in 



Fig. 1141 further elucidates the difference between the AMR 
and the static grid runs. Until t = 2 Gyr, the static run 
shows a slightly larger gas stripping than run A, similar to 
run D, indicating a more effective description of the onset 
of the KHI and of the following back-flow. In the late phase 
and especially after t = 3 Gyr this trend reverses, probably 
because the AMR runs fail to resolve the finest details of 
the dispersing structures in the subcluster tail. 

To sum up, the comparison with the static grid run 
shows that the new AMR simulations are able to correctly 
reproduce the main features of the subcluster evolution (see 
also Table [5} , though the fixed grid is more effective in the 
very late phase of the gas stripping. The benefit of AMR can 
be fully appreciated also from the computational point of 
view: in the case of idealised subcluster simulations, the use 
of such a static grid causes an additional computational cost 
of a factor from 15 to 40, depending on the AMR criteria, 
with increases also for the memory requirements and the 
output storage. 



6 DISCUSSION AND CONCLUSIONS 

We presented an interesting case of new refinement tech- 
niques, developed for resolving turbulent flows in AMR hy- 
drodynamical simulations, which may have useful applica- 
tions in astrophysical problems. The AMR is not only a 
tool used for this study, but has itself been investigated and 
tested in its capabilities. New reflnement criteria based on 
the regional variability of control variables of the flow were 
introduced and shown to be superior for simulations of tur- 
bulent media. 

As a flrst application, we presented a set of simulations 
of an idealised substructure in a wind in order to study a 
subcluster merger in a simplifled fashion. In this setup, the 
novel AMR criteria are suitable for the study of the tur- 
bulent wake, which is better resolved than using refinement 
criteria based on local gradients of selected hydrodynamical 
variables. The optimal control of the refining procedure en- 
sures high resolution where desired, as pointed out by the 
comparison with an equivalent static grid simulation. 

The overall evolution of the s ubcluster is similar t o 
the hydrodynamical simulations of lAcreman et all l|2003h . 
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Figure 13. Density slices of tiie subclustcr in tiie simulation with a 512'^ static grid. Left panel: t = 2 Gyr. Right panel: t = 3 Gyr. The 
density is linearly colour coded following the colour bar on the left. 



iTakizawal (1200531 ) and iTakizaTOl (|2005bl '). where the devel- 
opment of shear instability leads to the formation of a turbu- 
lent wake. The typical rms velocity in the wake is of the order 
of 500 km s^^, similar to th e above cited simula tions and to 
the theoretical predictions (|Subramanian et al .' 2006). This 
value is obtained only in runs using the new AMR criteria 
and tends to the result from an equivalent static grid simu- 
lation, whereas in the reference run A the inferred turbulent 
velocity is smaller by a factor of 2. 

We investigated the proper refining of the turbulent 
wake with several numerical experiments. The choice of the 
AMR criteria has an impact both on the turbulent velocity 
and on the morphology of the subcluster core as a result of 
the back-flow. We claim that this issue should be carefully 
considered in numerical simulations of subcluster mergers, 
as for instance in studies of cold fronts. If the problem is ad- 
dressed with an inadequate numerical setup (as in run A), 
the turbulent back-flow will be underestimated significantly. 
In extreme cases (run B, or Sect. l5.3TT|) . the simulation might 
fail to reproduce the flow downstream of the subcluster. We 
further showed that refining turbulence in the wake does not 
depend only on AMR but, as expected, also on the hydrody- 
namic solver. A detailed code compari son in this setup, lik e 
the one performed for the blob test bv lAgertz et al.l l|2007l ). 
would certainly be useful. 

As already stated, we do not intend to discuss the fea- 
tures of the cold front which forms ahead of the subcluster 
because our simulations do not inc lude the magnet i c field 
evolution and/or heat conduction (|Asai et al.l bOO'j . |2005| . 
I2OO7I : IXiang et al.|[2007l '). The observed cold fi-ont tempera- 
ture profiles (the best studied example is A3667) call for the 
suppression of the thermal conduction at the interface be- 
tween the cold subcluster gas a nd the hot ICM background 
l|Markevitch fc VikhhninlbOOTh . This suppression, in turn, 
would be naturally explai ned by a magnetic flux tan gential 
to the cold front surface. iDursi fc Pfrommeij l|2008l ) show, 
with MHD simulations using AMR, that the flow past a 
moving bullet generates vorticity, even if a magnetic fleld 
stabilises its surface against the KHI. The presence of MHD 
turbulence and further ampliflcation of the magnetic fleld 
is also suggested. Despite the simplifled approach of both 
setups, we infer that the generation of a turbulent wake and 
a proper modelling of this phase are therefore crucial also 



for the subcluster MHD simulations, and particular care in 
modelling the turbulent flow should be taken also in this 
class of simulations (with AMR or not). Only a clear as- 
sessment of the capability of a code in dealing with turbu- 
lent flows can allow to distinguish, for example, between the 
instability suppression caused by a magnetic field, and an 
unreliable modelling of the KHI onset. 

Our simplified approach to the minor merger case shares 
the limitations of previous, similar simulations cited above. 
However, it has the advantage of presenting a well controlled 
setup, which would be difflcult to study in detail in the 
framework of a full cosmological simulation. In view of the 
results it is necessary to note that, in the performed sim- 
ulations, the merging subcluster was embedded in a back- 
ground with a very simple velocity fleld. In a less idealised 
flow (cf. lNorman fc Brvan|[l999l : iDolag et 31.1120051 ). it is ad- 
mittedly unlikely that the turbulent tail would extend for 
several subcluster radii without being mixed. Nonetheless, 
it would not prevent from testing the turbulent charac- 
ter of the flow im mediately downstream of the subcluster 
(Takizawal l2005bl ). both observationally and in cosmological 
simulations. 

From an observational point of view, the turbulence 
level of the subcluster wake can be an interesting test for 
supporting the studied scenario. The rms velocity behind 
the subcluster can be used as an observable for turbulence 
for future X-ray spectrometers. This test would pertain to 
the wider, long- standing problem of detecting turbulence in 
galaxy clusters l|Schekochihin fc Cowlevll2006l and reference 
therein, for an overview). 

Although signiflcant efforts are un dertaken for mod- 
elling turbulent fl ows in SPH simulations (iDolag et al ] |2005l : 
IVazza et al.ll200^ ). grid-based AMR simulations are a very 
powerful approach to simulations in this field. The AMR cri- 
teria used in this work will be further applied to astrophys- 
ical problems, the next step being the simulation of galaxy 
cluster formation and evoluti on (c f . Pap er II) . According to 
the analysis of .Subramanian et al.l (|2006l ). besides the minor 
merger phase there are several other physical regimes when 
turbulence evolves and can play an important role in the en- 
ergy budget of a galaxy cluster. The present study corrobo- 
rates the importance of the minor mergers in the production 
of turbulence in galaxy clusters. Moreover, it provides new 
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tools for the theoretical study of this problem, which will be 
continued in Paper II. The final goal of our investigation is 
a consistent AMR modelling of tu rbulent flows based o n a 
subgrid scale model approach (cf. ISchmidt et af]|2006al lbh. 
This will be the topic of forthcoming work. 
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